function [data, crd, ftrf] = aft_SampleData3D(hfile, crd, opts)
% AFT::SampleData  - sample spatial data
%
% FORMAT:       data = obj.SampleData3D(crd [, opts]);
%
% Input fields:
%
%       crd         either a 4x3 range definition or Cx3 coordinates
%       opts        optional settings
%        .fmr       required for AMR/MAP objects
%        .mapvol    required for CMP/VMP objects, defaults to 1
%        .method    resampling method, one of {'cubic'}, 'linear', 'nearest'
%        .space     either of 'bvc', 'bvi', 'bvs', {'tal'}
%        .trans     4x4 transformation matrix, used for volume datatypes
%        .trf       TRF object(s) passed on to samplefmrspace
%        .v16       for VMR objects, if true access available V16 data
%
% Output fields:
%
%       data        either XxYxZ or Cx1 data
%
% TYPES: AMR, CMP, DMR, FMR, HDR, HEAD, MAP, MSK, VDW, VMP, VMR, VTC
%
% Note: for AMR/DMR/FMR/MAP objects, the trf option must be set with
%       at least the IA/FA (or FA alone) TRF object(s) to work

% Version:  v0.8a
% Build:    9102122
% Date:     Oct-21 2009, 10:11 PM CEST
% Author:   Jochen Weber, SCAN Unit, Columbia University, NYC, NY, USA
% URL/Info: http://wiki.brainvoyager.com/BVQXtools

% argument check
if nargin < 2 || ...
    numel(hfile) ~= 1 || ...
   ~isBVQXfile(hfile, true) || ...
   ~isa(crd, 'double') || ...
    isempty(crd) || ...
   (~isequal(size(crd), [4, 3]) && ...
    (ndims(crd) ~= 2 || ...
     size(crd, 2) ~= 3))
    error( ...
        'BVQXfile:BadArgument', ...
        'Invalid call to %s.', ...
        mfilename ...
    );
end
if nargin < 3 || ...
   ~isstruct(opts) || ...
    numel(opts) ~= 1
    opts = struct;
end
if ~isfield(opts, 'fmr') || ...
    numel(opts.fmr) ~= 1 || ...
   (~isBVQXfile(opts.fmr, 'dmr') && ...
    ~isBVQXfile(opts.fmr, 'fmr'))
    opts.fmr = [];
end
if ~isfield(opts, 'mapvol') || ...
    numel(opts.mapvol) ~= 1 || ...
   ~isa(opts.mapvol, 'double') || ...
    isinf(opts.mapvol) || ...
    isnan(opts.mapvol) || ...
    opts.mapvol < 1
    opts.mapvol = 1;
end
if ~isfield(opts, 'method') || ...
   ~ischar(opts.method) || ...
   ~any(strcmpi(opts.method, {'cubic', 'linear', 'nearest'}))
    opts.method = 'cubic';
else
    opts.method = lower(opts.method(:)');
end
if ~isfield(opts, 'space') || ...
   ~ischar(opts.space) || ...
   ~any(strcmpi(opts.space(:)', {'bvc', 'bvi', 'bvs', 'tal'}))
    opts.space = 'tal';
else
    opts.space = lower(opts.space(:)');
end
if ~isfield(opts, 'tal') || ...
    numel(opts.tal) ~= 1 || ...
   ~isBVQXfile(opts.tal, 'tal')
    opts.tal = {};
end
if ~isfield(opts, 'trans') || ...
   ~isa(opts.trans, 'double') || ...
   ~isequal(size(opts.trans), [4, 4]) || ...
    any(isinf(opts.trans(:)) | isnan(opts.trans(:))) || ...
    any(opts.trans(4, :) ~= [0, 0, 0, 1])
    opts.trans = [];
end
if ~isfield(opts, 'trf') || ...
   (~iscell(opts.trf) && ...
    (numel(opts.trf) ~= 1 || ...
     ~isBVQXfile(opts.trf, 'trf')))
    opts.trf = {};
elseif ~iscell(opts.trf)
    opts.trf = {opts.trf};
end
for tc = numel(opts.trf):-1:1
    if numel(opts.trf{tc}) ~= 1 || ...
       (~isBVQXfile(opts.trf{tc}, 'tal') && ...
        ~isBVQXfile(opts.trf{tc}, 'trf'))
        opts.trf(tc) = [];
    end
end
opts.trf = opts.trf(:)';
if ~isfield(opts, 'v16') || ...
   ~islogical(opts.v16) || ...
    numel(opts.v16) ~= 1
    opts.v16 = false;
end

% switch on filetype
sc = bvqxfile_getscont(hfile.L);
bc = sc.C;
ft = lower(sc.S.Extensions{1});

% for DMR/FMR/VTC objects, resolve transio, otherwise accept map
switch (ft)
    case {'amr'}
        if isempty(opts.fmr)
            error( ...
                'BVQXfile:MissingOption', ...
                'Option field .fmr must be set to access AMR in 3D space.' ...
            );
        end
        resetc = false;
        odata = zeros([size(bc.Slice(1).AMRData), numel(bc.Slice)]);
        for slc = 1:numel(bc.Slice)
            if istransio(bc.Slice(slc).AMRData)
                bc.Slice(slc).AMRData = resolve(bc.Slice(slc).AMRData);
                resetc = true;
            end
            odata(:, :, slc) = bc.Slice(slc).AMRData;
        end
        if resetc
            bvqxfile_setcont(hfile.L, bc);
        end
    case {'cmp'}
        nvol = numel(bc.Component);
        opts.mapvol = round(min(nvol, opts.mapvol));
        if istransio(bc.Component(opts.mapvol).CMPData)
            bc.Component(opts.mapvol).CMPData = ...
                resolve(bc.Component(opts.mapvol).CMPData);
            bvqxfile_setcont(hfile.L, bc);
        end
        odata = bc.Component(opts.mapvol).CMPData;
    case {'dmr'}
        if istransio(bc.DWIData)
            bc.DWIData = resolve(bc.DWIData);
            bvqxfile_setcont(hfile.L, bc);
        end
        switch bc.DataStorageFormat
            case {2}
                svol = 3;
            case {3}
                svol = 4;
            case {4}
                svol = 1;
            otherwise
                error( ...
                    'BVQXfile:IllegalSetting', ...
                    'Illegal DataStorageFormat setting in DMR object.' ...
                );
        end
        nvol = size(bc.DWIData, svol);
        opts.mapvol = round(min(nvol, opts.mapvol));
        switch (svol)
            case {1}
                odata = squeeze(bc.DWIData(opts.mapvol, :, :, :));
            case {3}
                odata = squeeze(bc.DWIData(:, :, opts.mapvol, :));
            case {4}
                odata = squeeze(bc.DWIData(:, :, :, opts.mapvol));
        end
    case {'fmr'}
        resetc = false;
        for slc = 1:numel(bc.Slice)
            if istransio(bc.Slice(slc).STCData)
                bc.Slice(slc).STCData = resolve(bc.Slice(slc).STCData);
                resetc = true;
            end
        end
        if resetc
            bvqxfile_setcont(hfile.L, bc);
        end
        szdt = size(bc.Slice(1).STCData);
        nvol = szdt(3);
        opts.mapvol = round(min(nvol, opts.mapvol));
        switch bc.DataStorageFormat
            case {1}
                odata = zeros([szdt(1:2), numel(bc.Slice)]);
                for slc = 1:numel(bc.Slice)
                    odata(:, :, slc) = bc.Slice(slc).STCData(:, :, opts.mapvol);
                end
            case {2}
                odata = squeeze(bc.Slice.STCData(:, :, opts.mapvol, :));
            otherwise
                error( ...
                    'BVQXfile:IllegalSetting', ...
                    'Illegal DataStorageFormat setting in DMR object.' ...
                );
        end
    case {'hdr'}
        odata = bc.VoxelData;
        if istransio(odata)
            odata = resolve(odata);
        end
        
        % also get transformation matrix and set space to bvc!
        cfr = hdr_CoordinateFrame(hfile);
        if isempty(opts.trans)
            opts.trans = inv(cfr.Trf);
        else
            opts.trans = opts.trans * inv(cfr.Trf);
        end
        opts.space = 'bvc';
    case {'head'}
        nvol = numel(bc.Brick);
        opts.mapvol = round(min(nvol, opts.mapvol));
        odata = bc.Brick(opts.mapvol).Data;
        if istransio(odata)
            odata = resolve(odata);
        end
        
        % also get transformation matrix and set space to bvc!
        cfr = head_CoordinateFrame(hfile);
        if isempty(opts.trans)
            opts.trans = inv(cfr.Trf);
        else
            opts.trans = opts.trans * inv(cfr.Trf);
        end
        opts.space = 'bvc';
    case {'map'}
        if isempty(opts.fmr)
            error( ...
                'BVQXfile:MissingOption', ...
                'Option field .fmr must be set to access MAP in 3D space.' ...
            );
        end
        resetc = false;
        odata = zeros([size(bc.Map(1).Data), numel(bc.Map)]);
        for slc = 1:numel(bc.Map)
            if istransio(bc.Map(slc).Data)
                bc.Map(slc).Data = resolve(bc.Map(slc).Data);
                resetc = true;
            end
            odata(:, :, slc) = bc.Map(slc).Data;
        end
        if resetc
            bvqxfile_setcont(hfile.L, bc);
        end
    case {'msk'}
        if istransio(bc.Mask)
            bc.Mask = resolve(bc.Mask);
            bvqxfile_setcont(hfile.L, bc);
        end
        odata = bc.Mask;
    case {'vdw'}
        if istransio(bc.VDWData)
            bc.VDWData = resolve(bc.VDWData);
            bvqxfile_setcont(hfile.L, bc);
        end
        nvol = size(bc.VDWData, 1);
        opts.mapvol = round(min(nvol, opts.mapvol));
        odata = squeeze(bc.VDWData(opts.mapvol, :, :, :));
    case {'vmp'}
        nvol = numel(bc.Map);
        opts.mapvol = round(min(nvol, opts.mapvol));
        if istransio(bc.Map(opts.mapvol).VMPData)
            bc.Map(opts.mapvol).VMPData = resolve(bc.Map(opts.mapvol).VMPData);
            bvqxfile_setcont(hfile.L, bc);
        end
        odata = bc.Map(opts.mapvol).VMPData;
        if bc.Map(opts.mapvol).EnableClusterCheck && ...
            isequal(size(odata), size(bc.Map(opts.mapvol).VMPDataCT))
            odata = odata .* bc.Map(opts.mapvol).VMPDataCT;
        end
    case {'vmr'}
        if ~isempty(bc.VMRData16) && ...
            opts.v16
            if istransio(bc.VMRData16)
                bc.VMRData16 = resolve(bc.VMRData16);
                bvqxfile_setcont(hfile.L, bc);
            end
            odata = bc.VMRData16;
        else
            if istransio(bc.VMRData)
                bc.VMRData = resolve(bc.VMRData);
                bvqxfile_setcont(hfile.L, bc);
            end
            odata = bc.VMRData;
        end
    case {'vtc'}
        if istransio(bc.VTCData)
            bc.VTCData = resolve(bc.VTCData);
            bvqxfile_setcont(hfile.L, bc);
        end
        nvol = size(bc.VTCData, 1);
        opts.mapvol = round(min(nvol, opts.mapvol));
        odata = squeeze(bc.VTCData(opts.mapvol, :, :, :));
    otherwise
        error( ...
            'BVQXfile:ObjectTypeUnsupported', ...
            'Unsupported object type for ::SampleData3D method.' ...
        );
end

% sample FMR based data differently
if any(strcmp(ft, {'amr', 'dmr', 'fmr', 'map'}))
    
    % take care of special coordinates
    if size(crd, 1) == 4 && ...
        all(isinf(crd(1, :)))
        crdx = cell(1, 3);
        crds = zeros(1, 3);
        crdx{1} = crd(2, 1):crd(3, 1):crd(4, 1);
        crdx{2} = crd(2, 2):crd(3, 2):crd(4, 2);
        crdx{3} = crd(2, 3):crd(3, 3):crd(4, 3);
        crds(1) = numel(crdx{1});
        crds(2) = numel(crdx{2});
        crds(3) = numel(crdx{3});
        [crdx{1:3}] = ndgrid(crdx{1}, crdx{2}, crdx{3});
        crd = [crdx{1}(:), crdx{2}(:), crdx{3}(:)];
    else
        crds = [size(crd, 1), 1];
    end
    
    % convert coordinates to BV system
    if strcmp(opts.space, 'tal')
        crd = 128;
    elseif strcmp(opts.space, 'bvi')
        crd = crd(:, [3, 1, 2]);
    elseif ~strcmp(opts.space, 'bvs')
        error( ...
            'BVQXfile:InvalidOption', ...
            'Sampling of DMR/FMR based data must not have bvc space.' ...
        );
    end
    
    % use AMR/MAP option DMR/FMR
    if ~isempty(opts.fmr)
        hfile = opts.fmr;
    end
    
    % sample data with samplefmrspace function
    ftrf = [opts.trf, opts.tal];
    data = reshape(samplefmrspace(odata, crd, hfile, ftrf, opts.method), crds);
    
% sample "already 3D" formats
else
    
    % get bounding box for valid objects and get transformation
    if ~strcmp(opts.space, 'bvc')
        bbox = aft_BoundingBox(hfile);
        ftrf = {bvcoordconv(zeros(0, 3), [opts.space '2bvc'], bbox)};
    else
        ftrf = {};
    end
    
    % any additional transformation?
    if ~isempty(opts.trans)
        if isempty(ftrf)
            ftrf = {opts.trans};
        else
            ftrf = {ftrf{1} * opts.trans};
        end
    end

    % manually apply transformation?
    if ~isempty(ftrf) && ...
       (size(crd, 1) ~= 4 || ...
       ~all(isinf(crd(1, :))))
        crd(:, 4) = 1;
        crd = (ftrf{1} * crd')';
        crd(:, 4) = [];
        ftrf = {};
    end

    % sample data
    data = flexinterpn_method(odata, crd, 0, ftrf{:}, opts.method);

end
